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Abstract 

O The accurate modeling of gravitational radiation is a key issue for gravitational wave astronomy. As 

T— I 

O simulation codes reach higher accuracy, systematic errors inherent in current numerical relativity wave- 

;-H extraction methods become evident, and may lead to a wrong astrophysical interpretation of the data. In 

^ this paper, we give a detailed description of the Cauchy-characteristic extraction technique applied to binary 

Q black hole inspiral and merger evolutions to obtain gravitational waveforms that are defined unambiguously, 

^ that is, at future null infinity. By this method we remove finite-radius approximations and the need to ex- 

O 

trapolate data from the near zone. Further, we demonstrate that the method is free of gauge effects and thus 
bJO is affected only by numerical error. Various consistency checks reveal that energy and angular momentum 

^ are conserved to high precision and agree very well with extrapolated data. In addition, we revisit the com- 

putation of the gravitational recoil and find that finite radius extrapolation very well approximates the result 
00 , 

^ at J'^. However, the (non-convergent) systematic differences to extrapolated data are of the same order of 

^ magnitude as the (convergent) discretisation error of the Cauchy evolution hence highlighting the need for 

^ correct wave-extraction. 

o 

> PACS numbers: 04.25 .dg, 04.30.Db, 04.20.Ha, 04.30.Nk 

X 
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I. INTRODUCTION 



The computation of gravitational radiation from black hole merger events has attracted con- 
siderable attention, since the pioneering work by Smarr and collaborators [[IH31. With the advent 
of ground-based laser interferometric gravitational wave detectors, as well as the prospect of the 
Laser Interferometer Space Antenna (LISA), interest in the problem has increased considerably in 
recent years. The measurement of gravitational waves is expected to provide an important probe 
of strong-field nonlinear gravity, thereby enabling empirical testing of this regime, as well as to 
provide direct observations of processes of fundamental interest in astrophysics. The sensitivities 
of LISA, and of the upcoming advanced ground-based detectors AdLIGO and AdVirgo, are high 
enough that even an error in aspects of the waveform calculation of 0.1% could lead to an incorrect 
interpretation of the astrophysical properties of a source, or of a test of general relativity flU. 

Nowadays, there are several codes that can produce a stable and convergent simulation of a 
black hole spacetime [|5]-fT4ll. and of particular topical relevance is the construction of long and 
accurate waveforms which can be used to analyze distinguishable features within the signals [fTSl , 
and also to construct a family of templates with the goal of improving gravitational wave detec- 
tion algorithms. Here the requirements are particularly challenging for numerical simulations, 
requiring waveforms which are accurate in phase and amplitude over many cycles to allow for 
an unambiguous matching to post-Newtonian waveforms at large separation. Some recent studies 
have shown very promising results in this direction for particular binary black hole models ffT6] - 

m. 

However, a particular difficulty arises from the fact that gravitational energy cannot be defined 
locally in general relativity, and is well-defined only at future null infinity, J'^. Since standard 
numerical evolutions are carried out on finite domains, there is a systematic error caused by esti- 
mating the gravitational radiation from fields on a world-tube at finite radius and the uncertainty 
in how it relates to measurement at J'^ [29|. Even though this error may be small, the expected 
sensitivity of upcoming detectors means that it is important to obtain an accurate result. 

A rigorous formalism for the global measurement of gravitational energy at J'^ has been in 
place since the pioneering work of Bondi, Penrose and collaborators in the 1960s [|30l l3T1l : and 
subsequently, techniques for calculating gravitational radiation at J'^ have been developed. The 
idea which we pursue here is to combine a standard Cauchy or "3 -i- 1" numerical relativity code 
with a characteristic code [[32ll in order to transport the wave information to J'^ using the full 
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Einstein equations. This method involves a Cauchy computation within a domain bounded by 
a world-tube, T, together with a characteristic computation with inner boundary T. If the data 
is consistently passed back and forth between the evolution schemes at T during the evolution, 
the method is known as Cauchy-characteristic matching (CCM). Given astrophysical initial data, 
such a method has only discretisation error [33], and a complete mathematical specification has 
been developed [|34|| . However, efforts to implement CCM encountered stability problems at the 
interface, and in the general case have not been successful. An alternative approach is Cauchy- 
characteristic extraction (CCE), or characteristic extraction. In this case, an independent Cauchy 
computation is carried out initially, with its usual timelike artificial outer boundary and which 
knows nothing about a characteristic code. Data on the world-tube F is stored, and then used to 
provide inner boundary data for a characteristic code - but not vice- versa. Initial implementations 
of CCE [|35ll36ll have shown the effectiveness of the method for a number of test problems. 

In this paper we describe the implementation of a CCE code, and present results obtained from 
the code for the real astrophysical problem of the inspiral and merger waveform of two black holes. 
The gravitational waveforms are calculated at , and are thus, for this case, the first waveforms 
that are unambiguous in the sense of being free of gauge and finite-radius effects. Further, the 
code is general purpose, in that it is independent of the details of the Cauchy code: it runs as 
a post-processing operation, taking as input only the required geometrical data on a world-tube. 
Thus it can be used in combination with other Cauchy codes, and applied to other astrophysical 
problems. 



The next section summarizes notation and results needed from other work. Then, Sec. Ill 



describes in some detail our implementation of characteristic extraction. Sec. IV presents results 
of testing the code against analytic solutions. The computation, including results, of the inspiral 
and merger of two black holes is described in Sec.|Vj Some technical details of the characteristic 



and extraction code parameters are described in Sec. VI 



n. REVIEW OF RESULTS NEEDED FROM OTHER WORK 
A. The Bondi-Sachs metric 

The formalism for the numerical evolution of Einstein's equations, in null cone coordinates, is 
well known Il30ll33l[37] - l40l . For the sake of completeness, we give a summary of those aspects of 
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the formalism that will be used here. We start with coordinates based upon a family of outgoing 
null hypersurfaces. We let u label these hypersurfaces, {A = 2,3), label the null rays, and 
r is a surface area coordinate. In the resulting = (m, r, x^) coordinates, the metric takes the 
Bondi-Sachs form [l30ll4ni 

-2e^^dudr - 2r'^hABU^dudx'^ + r'^hAsdx^dx'^ , (1) 

where h^^hsc = ^'^^ det{hAB) = (iet{qAB), with qab a metric representing a unit 2-sphere 
embedded in flat Euclidean 3-space; Wc is a normalized variable used in the code, related to the 
usual Bondi-Sachs variable by \^ = r + Wcv'^. As discussed in more detail below, we represent 
(lAB by means of a complex dyad qA- Then, for an arbitrary Bondi-Sachs metric, Hab can be 
represented by its dyad component 

J = hABq^q^'l^, (2) 
with the spherically symmetric case characterized by J = 0. We also introduce the fields 

K = + J J, U = U^qA, (3) 

as well as the (complex differential) eth operators 9 and B [|42l . 

The 10 Einstein equations Ra/3 = ^^{Tap — \gapT) are classified as: (i) hypersurface equa- 
tions - q'^RiA, h^^RAB - forming a hierarchical set for (3, U and W^, (ii) evolution equation 
q'^q^RAB for J; and (iii) constraints Roa- An evolution problem is normally formulated in the 
region of spacetime between a timelike or null world-tube, F, and future null infinity (J^^), with 
(free) initial data, J, given on n = 0, and with boundary data for j3, U, J satisfying the con- 
straints given on the inner world-tube. The inclusion of J7+ in the computational grid is made 
possible by compactifying the radial coordinate, r, for instance by means of a transformation of 
the form 

r 

r — )■ X = . (4) 

r + r^t 

In characteristic coordinates (but not in general), the Einstein equations remain regular at ^7+ 
under such a transformation. 
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B. The spin-weighted formalism and the S operator 



A complex dyad qa is a 2-vector whose real and imaginary parts are unit vectors that are 
orthogonal to each other. Further, qa represents the metric, and has the properties 

q^QA = 0, q-^QA = 2, qab = ^iqAQB + qAqs)- (5) 

Note that qa is not unique, up to a unitary factor: if qa represents a given 2-metric, then so does 
q'A — ^^""qA- Thus, considerations of simplicity are used in deciding the precise form of dyad to 
represent a particular 2-metric. The dyad commonly used to represent the unit sphere metric in the 
stereographic coordinates used here, is 



+ 1 + + 



Having defined a dyad, we may construct complex quantities representing all manner of ten- 
sorial objects, for example Xi = TAq^, X2 = T^^qAqs, Xs = Tq^ qAquff ■ Each object has 
no free indices, and has associated with it a spin-weight s defined as the number of q factors less 
the number of q factors in its definition. For example, s{Xi) = 1, s{X2) = 0, s^X-^) = —3, and, 
in general, s{X) = —s{X). We define derivative operators 9 and 5 acting on a quantity V with 
spin-weight s 

W = q^dAV + sTV, W = q^dAV - sTV (7) 
where the spin-weights of dV and dV are s + 1 and s — 1, respectively, and where 

r = -Iq^f^Aqs. (8) 

In the case of stereographic coordinates, T = q + ip. 

The spin-weights of the quantities used in the Bondi-Sachs metric are 

s{Wc) = siP) = 0, s(J) = 2, s{J) = -2, siK)=0, s{U) = 1, s(f7) = -1. (9) 

We will be using spin-weighted spherical harmonics [|43l l44ll ^Y^m (the suffix ^ denotes the 
spin-weight), using the formalism described in [45], and note that this reference gives explicit 
formulas for the sYim in stereographic coordinates. In the case s = 0, the s will be omitted, i.e. 
= oYim- Note that the effect of the S operator acting on Yim is 

gy,„ = ^i{e + i) d'Y^^ = ^{i-i)i{i + i){i + 2) ^y,^. (lo) 
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C. Gravitational radiation: News and V'4 



The mathematical theory relating metric quantities at future null infinity, , to gravitational 
radiation, and using the present formalism, is given in [|35ll38ll46l . In the original work of Bondi 
et al. [l30l . it was possible to assume that the coordinates at were such that /3 = J = there, 
and in that case the gravitational news takes the very simple form 

N=UudpJ (11) 

(where we have written the radial derivative in terms of p := rather than the usual notation, 
£, to avoid confusion with the spherical harmonic index t). However, the coordinates used in the 
characteristic code are fixed at the inner boundary F, and in general the Bondi gauge conditions 
are not satisfied. Previous work has presented the formalism for calculating, in a general gauge, 
the gravitational news [|38l . as well as Il35l . 

Geometrically, the Bondi gauge condition J = means that, for large r, a 2-surface of constant 
(m, r) is spherical rather than being, for example, an ellipsoid; and this condition can be expressed 
algebraically by saying that the 2-surface has constant curvature. An expression for the news, in a 
general gauge, must take (implicit) account of the transformation to Bondi gauge coordinates, 

r r^B]= ^{u^x'^y, xj^] = (12) 

The constant curvature condition leads to the factor u satisfying 

2K - mK + \{^^J + ^^J) + ^ {^J^J - ^J^J) =2{u^ + H^^DaDb log ou) (13) 

with all metric quantities evaluated at J^. Alternatively, we can initialize J = at m = 0, so that 
w = 1 at M = 0, and then evolve to by 

duUJ = -^ {Udu + UBu) -^{W + W) . (14) 

The Bondi gauge condition /3 = at ^7+ means that, for large r, coordinate time is the same as 
proper time, and the implementation of this condition is straightforward. 

While the transformation to Bondi gauge coordinates can be done explicitly B6l . in the present 
code it is done implicitly, leading to an expression for the news in terms of the code metric variables 
and coordinates. The formula is long and complicated [[38l and is not repeated here. However, 
in the linearized case, the news function is much simpler PTl . The results presented later are 
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approximately in the linearized regime, and so it can be instructive to examine this case (even 
though the code is fully nonlinear). The linearized formula for the news is 

N=U^dpJ+U\u + 2P). (15) 

In this case uj is simply related to J at , if the metric is decomposed into spin- weighted spherical 
harmonics, 

J = Jim 2^m, /3 = l^lm yim, (16) 

im Ira 

Then the news is 

tm (.ra 

(17) 

In characteristic work, it is conventional to work with a quantity \&4 rather than the usual V'4 
used in current Cauchy codes (see e.g.\9\). The relationship is 

Vl/4 = -(l/2)7/>4, (18) 

and we will usually transform the characteristic \I'4 to the standard ^4. The formulas for \I'4 are 
long and complicated and are not given here, but in the linearized case, \E'4 can be expressed in 
a simple way in terms of metric quantities at [|35ll (in which the notation ^ rather than \l/4 is 
used) 

^4 = \d%j - \d^j - \w - ^-^\m + m) + 9„g2/3. (i9) 

It can be shown that the linearized Einstein equations lead to the relation \E'4 = duN . 



D, Characteristic extraction 

We refer to [|34ll (see also [l36ll48ll ) for a full description of the characteristic extraction process. 
While characteristic extraction appears to be only a coordinate transformation, it is actually rather 
more complicated. The difficulty is that Bondi-Sachs coordinates use a surface area radial coordi- 
nate, and this coordinate can be constructed only once the angular coefficients of the metric have 
been found. Thus the construction proceeds in two stages. 

In the first stage, we use an affine coordinate A in the radial direction, and find the transformed 
metric and its first a derivatives at angular grid-points of the extraction world-tube V. The process 
is illustrated in Fig.[T] and can be summarized as follows: 
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FIG. 1: Schematic illustration of the (first stage) construction of characteristic coordinates and metric. 

• Define a world-tube Thy x'^+y'^+z'^ — R"^ with R constant, and induce angular coordinates 
0^ on r as though in Euclidean space. 

• Let H be a hypersurface of constant t, and define the 2-surface S — H (IT. 

• Let be a unit normal to H, and let be normal to S in H. 

• Construct outgoing null vectors = + s", and then outgoing null geodesies in direction 

with afflne parameter A. 

• The union of such geodesies is the null cone N, labeled by u being the Cauchy time t where 

meets T. 

• Construct the Jacobian for the coordinate transformation {t, x, y, z) — )■ {u, A, (j)^). 

• Find the transformed metric at the angular grid-points of F. 

• Find the first ,a derivatives of the transformed metric at the angular grid-points of F. 

In the second stage, we make the transformation to a surface area coordinate r. The difficulty 
here is that, in general, r is not constant on F. Thus in order to set data on an inner world-tube 
of the characteristic grid, we need the metric quantities on the world-tube, as well as their first 
derivatives off the world-tube. The process can be summarized as follows: 

• Make the coordinate transformation A ^ r = r{u, A, 0^) with r defined by the condition 
that it is a surface area coordinate. 
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• Find r and r a at the grid points of T. 

• Find the metric and its first a derivatives at the grid points of T. 

• Find J,U,Wc and j3 at the grid points of F. 

• Using dr = d\/r^x, find J.,., U^r, W^c,r and [5,^ at the grid points of F. 

• Set J, U, U^r, Wc and (3 on an inner world-tube of the characteristic grid. 

The value used for f/^. is second-order accurate on the world-tube F, and it is this value which is 
used at the inner world-tube of the characteristic grid, so, in general, it is only first-order accurate. 
However, in practice, the overall performance of the code remains second-order accurate. 

E. The Cauchy evolution code and finite radius measurements 

A fundamental problem with early attempts to implement CCM was the sensitivity of the 
3-1-1 (Cauchy) implementation to boundary perturbations. In recent years, however, a number 
of strongly hyperbolic formulations have been developed and proven to be numerically robust 
for binary black hole evolutions using artificial outer boundary conditions We evolve the space- 
time in the neighborhood of the binary black holes using a variant of the "BSSNOK" evolution 
system Il49] - l5ll . according to the implementation described in [jH. This formulation uses a con- 
formally transformed metric, ^yab, and (traceless) extrinsic curvature, Aab, as evolution variables, 
where the conformal transformation is defined by 

= ^ In det 7a6, %b ■= e~^'*'-fab- (20) 
Instead of evolving cf) directly, we use the concomitant 

06:=(det7a6)"'/' = e-2^ (21) 

following the suggestion of [l52l . 

Initial data for binary black holes is determined using the puncture method [53], and corre- 
sponds to a conformally flat solution of the Hamiltonian constraint with the extrinsic curvature 
specified according to Bowen and York [1541 . Parameters for quasi-circular orbits starting from a 
given separation are determined by a post- Newtonian estimate [[T4]| . resulting in a low eccentricity 
(e ~ 0.04) inspiral for the equal-mass non-spinning binary considered below. 

9 



The lapse and shift are evolved according to the "1+log" [l55l and "f -driver" [|56l . respectively, 
including shift- ad vection terms which have proven to be important for maintaining phase accuracy 
in the puncture motion [I71I571- This combination of gauges are key ingredients of the "moving 
puncture" method and the control of evolution of the coordinate singularity within the black 
hole throat llSSll . 

The fields are discretize on regular locally Cartesian grids which cover the neighborhood of 
the binary to some finite coordinate radius. Finite differences at 8th-order are used to numerically 
approximate derivatives. Boundaries between coordinate patches are determined by 5th-order La- 
grange polynomial interpolation. Fields are evolved in time according to the BSSNOK prescrip- 
tion, using a 4th-order Runge-Kutta integrator. A detailed description of the numerical scheme 
can be found in flU. 

An objective of this paper is to compare the gravitational waveforms measured at J7+ with 
the locally determined measures which are in common use in numerical relativity. These local 
measurements generally evaluate the geometrical fields (the metric and associated derivatives) on 
a coordinate sphere of constant radius. The measurement spheres are assumed to be in the "wave- 
zone", which is commonly taken to be from r = 18M ll59l to r = 225M |l8l from the source. The 
particular numerical grid structure employed in this paper follows that of companion studies in 
which we have demonstrated accurate and convergent wave measurements to r = lOOOM [|9ll60l. 
Polynomial extrapolations to spatial r — t- oo are applied in order to reduce the systematic error 
associated with the finite radius measurements. Below, we demonstrate that the error associated 
with this procedure in comparison with the CCE result is, however, still larger than the numerical 
error for the models we have studied. 

Two methods are used in order to estimate the gravitational radiation at a given radius. The first 
evaluates the Weyl tensor in a particular null frame (n, 1, m, m) [31], which is chosen to be ori- 
ented along the coordinate radial direction, and orthonormalised using the local metric. According 
to the peeling property of asymptotically flat spacetimes, the component 

1p4 .= Cnmnm (22) 

is associated with the outgoing gravitational radiation. As an alternative to tp^, we also measure 
waves by evaluating the Zerilli-Moncrief variables, assuming the spacetime metric at large radius 
from the source to be representable as perturbations of a fixed Schwarzschild background [[6TI - 
[63l . We evaluate Ist-order gauge invariant odd-parity and even parity multipoles, (Q^^ and Ql^, 
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respectively) [|64l [65l . The implementation of the ip4 and Zerilli-Moncrief measurements follow 
the outline described in [|9l[T0l. 

Our evolution code is built using the Cactus computational framework [|66l |67l . The grid 
structures are implemented using extensions to the Carpet mesh refinement driver Il68l470ll . al- 
lowing the evolution domain to be covered by multiple patches with generalized local coordi- 
nates [9]. Initial data is computed using the TwoPunctures code [71]), which have been used 
for previous Cartesian evolutions. 

III. IMPLEMENTATION OF CHARACTERISTIC EXTRACTION AS POST-PROCESSING 
CODE 

The characteristic extraction is independently performed as a post-processing step after the full 
Cauchy evolution. During the Cauchy evolution, we output the Cauchy 4-metric and its deriva- 
tives. The variables are decomposed in terms of spherical harmonic amplitudes on a set of extrac- 
tion spheres with fixed radii around the origin of the black hole spacetime. On the characteristic 
side, this data can then be read and reconstructed to generate inner boundary data on the world- 
tube. 

Since we obtain the fields at ^7+ in terms of coordinate time u, we finally have to transform the 
fields to constant Bondi time ub- Furthermore, due to the rather complicated procedure of calcu- 
lating the conformal factor to for the calculation of the news, we have implemented a linearized 
approximation that turns out to be numerically more accurate than the full non-linear computation. 

In the next subsections, we will describe these procedures in more detail. 

A. World-tube boundary data 

Boundary data for the characteristic evolution are constructed during Cauchy evolution in terms 
of the ADM variables on coordinate spheres with fixed radius R. The time series of the coordinate 
spheres at successive times generates the discretize world-tube. 

In order to be able to transform the ADM variables to a representation of the full 4-metric in 
terms of Bondi coordinates and in terms of the characteristic evolution variables in a neighborhood 
of the world-tube, we have to compute the following quantities: the ADM metric components 7jj, 
as well as the time-derivatives dt 'jij and Cartesian derivatives dk 'Jij, the lapse a, as well as time 
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and Cartesian derivatives dt a and di a, the shift vector l3\ and again time and Cartesian derivatives 
dt /3* and di (3'^. In the specific formulation of the evolution system that is used in this paper, we 
can compute these quantities as follows. The BSSNOK evolution variables are transformed to the 
ADM variables after each evolution step and are therefore known everywhere. What remains is 
the calculation of time and spatial derivatives. 

The time derivatives of the lapse and shift vector are known directly from the RHS of the 
"1+log" slicing condition and the hyperbolic F-driver condition (see @| for details of the imple- 
mentation of these gauge conditions). The RHSs of these equations have to be computed at each 
time-step in order to evolve the gauge variables with a method of lines (MoL) scheme, similar to 
the evolution equations. 

The time derivatives of the metric components can be obtained via the ADM relation 

dt lij = -2aKij + + Djl3i , (23) 

where Di is the covariant derivative operator compatible with the spatial 3-metric 7ij and Kij 
the extrinsic curvature. In order to compute the covariant derivative, it is necessary to invert the 
3-metric and to calculate the Christoffel symbols on each grid point. Because the shift vector 
is given in contravariant form, we have to transform it to its covariant representation with index 
down. After looping over all points on the computational grid, we store the outcome in grid 
functions, i.e., one per metric component. 

Next, obtain the Cartesian derivatives of the metric, lapse and shift. This is straight-forwardly 
done by applying finite difference operators to the variables. Note that in case of non-Cartesian 
local coordinates, the derivatives are calculated via the global derivatives as described in dH. For 
practical reasons (particularly storage and memory requirements), it is sufficient to calculate the 
radial instead of the Cartesian derivatives. The Cartesian derivatives can later be reconstructed in 
terms of angular derivatives of the spherical harmonics which are known analytically. This saves 
two additional grid functions per ADM variable. We obtain the radial derivatives in terms of the 
Cartesian derivatives via 

dr = ^{xd,+ydy+zd,). (24) 

This completes the first step of calculating all variables that are necessary for the construction of 
characteristic boundary data. What remains is the projection of these variables on to a coordinate 
sphere with fixed radius R and a subsequent decomposition in terms of scalar spherical harmonics. 
To project the variables on to the sphere, we use fourth-order Lagrange interpolation. Afterwards, 
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the variables are decomposed as 

U^= [ dnY,^f{n), (25) 

for all variables / whether scalar, vector or tensor quantities. The projection and decomposition 
are both done on polar spherical coordinate spheres where the surface integration is done with a 
Gauss quadrature scheme as described in [9], and the resulting array of spherical harmonic modes 
is stored in a file for further processing. The reconstruction from the data in that file will be 
described in the next subsection. 



B. Reconstruction from harmonic modes 



Once the boundary data in terms of the ADM variables and their time and radial derivatives 
are known on a coordinate sphere with radius R and stored in a file, the characteristic boundary 
module reads and interprets the data in that file in order to reconstruct the variables on S"^. That is, 
the characteristic code is run as a post-processing tool to obtain the gravitational radiation at J'^. 
The variables are reconstructed via 

/ = ^ ftrrXtm, ■ (26) 

for any variable /. 

The advantage of reconstructing the variables from harmonic modes is the independence in 
terms of angular coordinates,' as well as resolution. Additionally, the necessary finite ^-mode cut- 
off can act as as a filter that factors out angular high-frequency noise. In practice for the binary 
black hole models considered here, the amplitudes of higher £ modes falls off rapidly so that they 
can be considered essentially irrelevant numerically for, say, £ > 10. 

Note that for practical reasons, we have preferred to calculate the radial instead of the Cartesian 
derivatives of the ADM variables. As a consequence, we have the harmonic modes of the radial 
derivatives. In order to obtain the Cartesian derivatives instead, one can take angular derivatives of 
the spherical harmonics and then apply a Jacobian transformation from stereographic to Cartesian 
coordinates 

if = ^ Or flmyiMP) + femTT-^ + Jlm^^ , (27) 



X* ox^ ox^ oq ox^ op 



^ On the Cauchy side, the ADM quantities are represented in polar- spherical coordinates on S^. The characteristic 
code, on the other hand, works with a stereographic coordinate mapping of S^. 
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where g,p denote stereographic angular coordinates. By using the relation 



_dA + QA- 2ipsA _ i{QA -QA + 2qsA) 

1 + + 1 + + P'^ 



for a quantity A with spin- weight s, we can express partial derivatives with respect to angular 
coordinates in terms of the eth-derivative 

d dr 

-f = TT— dr flmXlm 



dq /5Ygm + §>£mA , , dp . (Wim-^Ye,m\ 
+ Jem^— I — ^— — K- + jem^—i I -r- — ^— — ^ , K^y) 
ax* \ 1 + + p"' / ax* \ 1 + / 

where we have used the fact that the scalar spherical harmonics have spin-weight s = 0. This can 
be further simplified by replacing the eth-derivatives of the spherical harmonics by spin-weight 
s = ±1 spherical harmonics. Finally, we can write 



(9 f ^™ A /£(£ + l^l 

= 2^ 2 l^^^rn{q,i ' lV,i) ' -lYimil, + m)] + r.^Yim , (30) 

a X* 1 + g"' + p"' 



where 




(31) 

(r ± 2; — x^/r, —xy/r, —xz/r =p x) , (32) 



(r ± zY 

P,i = di I — ^ I = (^^^^y {Txy, ±r + zT y^/r, -y T yz/r) , (33) 

where the upper sign is valid for the north patch and the lower sign is valid for the south patch. 
This means that we do not have to calculate any extra derivatives numerically and by virtue of 



Eq. (30), we immediately obtain the Cartesian derivatives from the radial derivatives. 



Once the ADM variables and their derivatives are reconstructed on the coordinate sphere, the 



code executes the computations discussed in Sec. II D to obtain the boundary data in terms of the 
characteristic evolution variables and in terms of the Bondi coordinates. Once this step is complete, 
the code can use the boundary data to evolve the characteristic evolution variables further in time. 

Results for the application of this method to a binary black hole merger are presented in SecfV] 
However, the fields at first have to be transformed to constant Bondi time ub- This is described 
in the next subsection. 
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C. Interpolation to constant Bondi time and mode decomposition 



After the evolution of the characteristic variables to , the news function and the Weyl 
scalar ^4 are computed in terms of inertial spatial Bondi coordinates {rB,PB, Qb), but are not yet 
given in terms of the inertial Bondi time coordinate ub, i.e., we have N = N{rB, Pb, Qb,u). Also, 
one deals with the complete 2D data at but for analysis purposes, it is more convenient to 
have the data decomposed in terms of spherical harmonic modes. 

In a previous implementation, the transformation to constant ub was achieved by storing the 
news function as a function on S*^ at each time-step to a file. In a post-processing step, these time- 
slices were then read and interpolated to constant Bondi time. This, however, is rather inconvenient 
to use in practice. Rather, one would like to have the variables automatically transformed to the 
correct time coordinate and decomposed in terms of spherical harmonic modes, avoiding the need 
for additional post-processing tools. For this reason, we have extended the module for calculating 
the news function and the Weyl scalar \t'4 to take care of these issues. 

During characteristic evolution, we know the inertial Bondi time in terms of coordinate time 
and angular coordinates ub = UB{u,p,q) at each time-step at J'^. We have found that in practice 
the difference between u and is very small, so that by keeping five time levels of ub, N and ^4 
in memory at each time-step, it is possible to interpolate and 4 to ub = const, for each point 
on S"^. This is done by means of fourth-order Lagrange interpolation on the fly during evolution. 
In practice, we average the Bondi time over S*^ on the past-past time-level to find the target time 
Ub = const. We use the past-past time-level to have the interpolation stencil centered around the 
target interpolation time to maintain maximal accuracy in the interpolation. Once the target time 
is known, and \l/4 are interpolated for each point on S*^ to Ms = const, so that we finally have 
N = N{rB,PB,qB,UB) and ^4 = ^ A{rB,PB,qB,UB). 

Afterwards and if requested by the user, the code can automatically decompose the two func- 
tions in terms of spin- weighted spherical harmonic modes. Therefore, the data is available in a 
format that is convenient for analysis without any further post-processing. A set of post-processing 
tools can read in the harmonic modes of \l/4 and calculate radiated energy, linear and angular mo- 
mentum as well as the wave-strain h. 

However, there is one word of caution required. As mentioned earlier, the characteristic code 



computes a quantity \[^4 which is related to the commonly used ipA, via Eq. ( 18 1. It follows that the 
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harmonic modes are related by 

^ira ^ -2¥f'^{-l)'^ . (34) 

D. The linearized conformal factor 

In order to compute the news and ^4, it is first necessary to calculate the conformal factor 
oj of the conformal compactification of the Bondi metric at . However, the full non-linear 
computation of this factor can display a high level of noise when it is implemented numerically, 
reducing its accuracy. Since in all cases that we consider the fields at are in a linear regime, 
we can make use of a relation between the conformal factor u and J at in the linearized 
approximation. 

Given J at in terms of s = 2 spin- weighted spherical harmonics, i.e., 

J|j7+(m,X^)= ^ Jl,m{u)2yt,m, (35) 
l>2,m 

we can write uj in terms of scalar spherical harmonics as 

a;(M,X^) = 1 + X] iOi^rn{u)Yi^m, (36) 
l>2,m 

where 

The linearized computation has been implemented as an alternative to the full non-linear com- 
putation of oj and either method can be selected at run-time. 

IV. TESTS OF THE CHARACTERISTIC EXTRACTION CODE 

In this section, we report on tests we have performed in order to calibrate the characteristic 
extraction code. For this, we put an analytical solution to the ADM metric variables on the Cauchy 
side, generate the associated spherical harmonic modes on each timestep and afterward run the 
characteristic extraction code. 

A. Linearized solutions 

We have tested convergence of the entire extraction process against an analytical (linearized) 
solution [|47l by transforming this solution to the ADM variables, i.e., lapse a, shift /3' and the 3- 
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metric components 7ij, on the Cauchy side. According to the description in Sec. Ill we decompose 
the ADM variables on a world-tube after each timestep in order to produce spherical harmonic 
modes that can be used as boundary data for the subsequent characteristic evolution. On the 
characteristic side, we measure the error, that is, the difference between evolved variables and 
exact fields. We do so by performing a set of two resolutions, and we can assume convergence 
against the linearized solution given that the amplitudes of the fields are small, i.e., they are in the 
hnearized regime. 

The solutions to be used are those given in Sec. 4.3 of ref. PTll for the case of a dynamic 
spacetime on a Minkowski background with i = 2, m = 0. We write 

J = ^{t - IW +m + 2) 2F,^3?( J,(r)e^^"), 

= Y,^n{W^[r)e'-^), (38) 

where Ji{r\ Ui{r), (3^, Wce{r) are in general complex, and taking the real part leads to cos(z/m) 
and sm(i'u) terms. The quantities (3 and Wc are real; while J and U are complex due to the terms 
2Yim and lYem, representing different terms in the angular part of the metric. In the case i = 2 

(32 = /3o, 

24/3o + SiuCi - %v^C2 Ci C2 
-^^^'^ - 36 ^Tt~V2?>' 

. , -242z//3o + 3z/^Ci - v^C2 2/3o Ci ivC^ C2 
^2l^J 36 + ^ + 2r2 + 3r3 ^ Ar^' 

^-W = ^ + ^ + -r + ^' (39) 

with the (complex) constants (3q, Ci and C2 freely specifiable. 
We present results for the case £ = 2, m = 0, z/ = l, with 

Ci = lo-^ C2 = 3 X lo-^ /3o = o (40) 



inEq. (39). 



The numerical simulations use a compactified radial coordinate, Eq. Q, with r„t = 9. Data is 
prescribed at time u = and extracted at a number of different world-tube locations, specifically 
at r = {30, 60, 90, 120}. The stereographic grids (with ghost zones excluded) are 

Coarse: = Uq = Up = 41, Fine: = rig = Up = 81; (41) 
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FIG. 2: First-order convergence in the L2-norm of J at against a linearized solution as given by (39 1 



and there is no overlap between the two patches, i.e., we set the code parameter ggj^e ~ ^ which 
means that on the nominal grid, the holomorphic coordinate function = q + ip takes values in 
q,p E [—1, 1]. The compactified coordinate takes values x G [xin, 1] such that the world-tube is 
located on the characteristic grid and close to the inner boundary to avoid accuracy losses. For 
example, for a world-tube located at r = 30M, we set xin = 0.74. For different world-tube 
locations, Xin may be changed accordingly. 

The grids on the Cauchy side have a radial resolution of h = O.AM and h = 0.2M with the 
number of angular points per direction and per patch set to n = 21 and n = 41, respectively. The 
timestep size is taken to be At = O.IM and At = 0.05M. 

In Fig.|2} we show convergence in the L2-norm of J at J'^ to the linearized solution, starting off 
from a world-tube at r = 30M. We observe clear first-order convergence. This is consistent with 
the fact that the boundary data for U^r is only known to first-order off the world-tube. However, 
for other models presented in subsequent sections, it appears that the coefficient of numerical 
error arising from f/ ^ is small enough that a second-order convergence exponent can be measured, 
corresponding to the interior finite difference accuracy of the characteristic code. 



B. Gauge invariance of a shifted Schwarzschild black hole 

An important aspect of CCE as a wave extraction method is the fact that it should produce 
invariant results independent of the gauge at the world tube. As a simple test problem which 
demonstrates this property, we place a single static Schwarzschild black hole at the origin of our 
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spacetime, and apply a shift to the solution, according to 

= Au cos(wt) , = = , (42) 

so that the x-coordinate oscillates as 

x = A sin(a;t) , (43) 

where the amplitude is set to A = 1.0 and the frequency is set to w = 0.1. As the spacetime 
is spherically symmetric and static, the resulting radiation content must be zero, and hence, in 
numerical simulations the residual in and the news N must converge to zero. This is indeed 
what we find. Fig. |3] shows the {i, m) = (2, 0) mode of \E'4 and at two resolutions h = AM 
and h = 0.2M scaled for second-order convergence. All other modes converge to zero at the same 
order. 

In contrast to this, finite radius extraction methods based on tp^ and on gauge-invariant per- 
turbations of Schwarzschild fail to be convergent. This is shown in Fig. |4| where we plot the 
perturbative Zerilli-Moncrief (even) master function Qj^2 m=o 1164117211731 and -ipA for the two res- 
olutions h = 0AM and h = 0.2M without any rescaling. The two waves are on top of each other 
indicating that the results do not converge to zero. 

Although the effect is exaggerated by this artificial test, the point remains that the commonly 
used finite radius extraction methods are susceptible to gauge variations which are difficult to dis- 
entangle from the genuine wave signal. The CCE method as implemented removes this systematic 
error completely. 

V. RESULTS FOR BINARY BLACK HOLE INSPIRALS 

We have carried out full three-dimensional, non-linear evolutions of two binary black hole 
systems for the whole spacetime out to using the Cauchy and characteristic formulations as 
described in the previous sections. The first configuration represents an equal-mass non-spinning 
binary, with modes up to £ = 8 of the CCE waveform illustrated in Fig. [5| while the second 
configuration is an equal-mass binary with individual spins that are aligned with the orbital angular 
momentum and hence the system is not subject to precession effects. In this section, we compare 
the differences in the waveform measures at finite radius and at J'^, and compute various features 
of the two binaries. 
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FIG. 3: Second-order convergence of the gravitational wave signal to zero for the shifted-Schwarzschild 
test case. The left panel shows the £ = 2, m = mode of the Weyl scalar ^4 for the two resolutions 
h = 0.40M and h = 0.20M where the latter is scaled for second-order convergence. The right panel 
shows the {i, m) = (2, 0) mode of the news A'^ for the same resolutions and the same rescaling. Similar 
plots for all other modes exhibit the same order of convergence. 




t/M tjM 



FIG. 4: Non-convergent behavior of finite-radius measures of the gravitational radiation for the shifted- 
Schwarzschild test case. Shown are the perturbative Zerilli-Moncrief (even) master function {left) 
and 3^V'4 computed at a finite radius {right) extracted at r = lOOM. For each case, we plot two resolutions 
h = 0.20M and h = 0.40M of the {£, m) = (2, 0) mode without any rescaling. We find that the waveforms 
do not converge to zero, as should be expected for this test case. 

We demonstrate that for these tests, the wave-extraction at contains only numerical error. 
Additionally, we measure the energy, linear and angular momenta radiated through gravitational 
waves, both at finite radii and at and test for conservation of energy and angular momentum. 
As the emission of gravitational radiation is symmetric in the first case, but anti- symmetric in the 
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FIG. 5: Dominant spherical harmonic modes of V'4 for £ < 8 extracted at for the model ai = 02 = 0. 

second case, the net-linear momentum carried by the gravitational waves is non-zero for the second 
binary. Hence, this binary gives rise to a gravitational recoil, or kick, of the merger remnant. We 
measure this effect to high precision using gravitational waveform data at finite radii and at . 

A. Initial data 

For the first equal-mass non- spinning binary, the numerical evolution starts from an initial 
separation djM = 11.0, and goes through approximately 8 orbits (a physical time of around 
1360M), merger and ring-down. The masses of the punctures are set to m = 0.4872 and are 
initially placed on the x-axis with momenta p = (=fO. 000728, ±0.0903, 0), giving the initial slice 
an ADM mass Madm = 0.99051968 ± 2 x 10^^. The z-axis is chosen to coincide with the initial 
orbital angular momentum and is given in Table |Ij The y-axis is then chosen perpendicular to the 
other axis. 
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The second binary has been investigated mainly so as to be able to measure the kick-velocity 
that the merger remnant acquires due to the asymmetric emission of gravitational radiation. As 
this effect depends largely on the very last orbit (compare e.g. |[TOl ). it is sufficient to consider 
a binary with a closer initial separation which in our case has been chosen to be d/M = 4.0. 
The initial masses of the punctures are set to m = 0.2999 and are initially placed on the x- 
axis with momenta p = (=fO. 00211428, ±0.112539, 0), giving the initial slice an ADM mass 
Madm = 0.98826246 ± 2 x 10^^. 

The coordinate system is chosen such that the initial orbital angular momentum is aligned with 
the z-axis of the coordinate system, and we represent the initial spin vectors by 

Si = (0,0,ai), S2 = (0,0,a2), (44) 

where we have chosen (01,02) = (0,0) for the first binary and (01,02) = (0.8, —0.8) for the 
second one. 

The initial data parameters were determined using a post-Newtonian evolution from large initial 
separation, following the procedure outlined in [fT4l . with the conservative part of the Hamiltonian 
accurate to 3PN, and radiation-reaction to 3.5PN. The measured eccentricity of the resulting orbits 
is 0.004 ± 0.0005 

The most convenient way to set the initial data on the initial nuU-hypersurface of the charac- 
teristic grid is to set J(u = 0) = everywhere, since the wave-extraction algorithm requires 
J{u = 0)\j+ = 0. Fortunately, it turns out that J = on the extraction world-tube F at m = is 
consistent with conformal flatness, although in general there will be an inconsistency at points on 
the null cone off the initial Cauchy hypersurface. 

B. Grid setup 

The binary black hole evolution was carried out on a 7-patch grid structure, as described in 
||9l|60l, incorporating a Cartesian mesh-refined region which covers the near-zone, and six radially 
oriented patches covering the wave-zone. The inner boundary of the radial grids was placed at 
rt = 35. 2M relative to the centre of the Cartesian grid. The outer boundary for the spherical grids 
was chosen based on the expected time duration of the measurement and radius of the furthest 
extraction world-tube, in order to remove any influence of the artificial outer boundary condition. 
In particular, given that the evolution takes a time for the entire inspiral, merger and ringdown, 
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FIG. 6: Portions of the Kruskal diagram which are determined by the numerical evolution. The green 
horizontal lines indicate the region of spacetime that is determined by the Cauchy evolution and which has 
finite spatial extent with artificial outer boundary at r^,. The red diagonal lines indicate the region that is 
determined by characteristic evolution, and which has an inner boundary at a world-tube F located at rp, 
consisting of data from the intersecting slices (green horizontal lines) of the Cauchy evolution. Note that 
due to the final spatial extent of the latter, the future Cauchy horizon of the Cauchy initial data, indicated 
by the dotted diagonal line parallel to , principally prohibits the determination of the spacetime out to 
timelike infinity i+. However, as long as the world-tube F is located within the future Cauchy horizon, 
the numerically evolved subset of the spacetime is determined consistently. Particularly, the artificial outer 
boundary is sufficiently removed such that the evolved binary black hole spacetime forms a merged black 
hole before the Cauchy horizon reaches the world-tube. 

and a world-tube location at a finite radius rr, we would like to ensure that a disturbance traveling 
at the speed of light c = 1 [174117511 from the outer boundary does not reach the world-tube radius 
(see Fig. [6])^. Thus we place the boundary at 

rfe > + 2rr. (45) 

For the particular evolutions considered here, for the non-spinning binary and for the spinning 
binary respectively, we have evolution times ~ 1350M and — 550M, and outermost 

^ The 1 + log slicing condition which we use propagates at \f2c iTSSl, however this is a gauge mode and empirically 
we find it to have negligible effect on finite radius measurements. Particularly, since CCE is gauge invariant, this 
does not affect the waveforms at . 
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extraction world-tubes rr = lOOOM and rr = 500M. We have placed the outer boundary of the 
evolution domains at rf, = 3600M and r?, = 2000Af , respectively. 

The near-zone grids incorporate 5 levels of 2:1 mesh refinement, covering regions centred 
around each of the black holes. For the highest resolution we have considered here, the finest grid 
(covering the BH horizon) has a grid spacing of dx = 0.02M. The wave-zone grids have an inner 
radial resolution which is commensurate with the coarse Cartesian grid resolution, dr = 0.64M 
in this case. This resolution is maintained essentially constant to the outermost measurement 
radius (r = lOOOM), at which point we apply a gradual decrease in resolution (see BH) over a 
distance of r = 500M. From r = 1500M to the outer boundary, we maintain a resolution of 
dx = 2.56M, sufficient to resolve the inspiral frequencies of the dominant (i,m) = (2, 2) mode 
of the gravitational wave signal. The angular coordinates have 3 1 points (30 cells) in p and a on 
each of the 6 patches. The time-step of the wave-zone grids is dt = 0.144, and we take wave 
measurements at each iteration. 

We have carried out evolutions at three resolutions, h = 0.64M, h = 0.80M and h = 0.96M, 
in order to estimate the convergence of our numerical methods (compare [9]). 

The characteristic extraction world-tubes F have been placed at rr = lOOM and rr = 200M 
for the first binary, while for the second binary we have chosen rr = lOOM and rr = 250M. We 
note that over the course of the evolution, these world-tubes are causally disconnected from the 
artificial outer boundary. 

C. Invariance with respect to world-tube location 

Invariance with respect to the world tube location is demonstrated in Fig. |7j We have con- 
sidered the differences between waveforms at J'^ resulting from two independent characteristic 
evolutions using boundary data at rr = lOOM and rr = 200M , respectively, and for two reso- 
lutions, h = 0.96M and h = 0.64M. The difference between the results should be entirely due 
to the discretisation error, and indeed this is what we find. The differences converge to zero with 
approximately second-order accuracy, as expected for the null evolution code. The figure displays 
the differences in the amplitude and phase of the wave mode 'ip'ii^ = 2, m = 2) for resolution 
h = 0.96M and h = 0.64M scaled for second-order convergence. The same order of convergence 
is also obtained for higher order modes such as ^4(£ = 4, rri = 4) and ipi{i = 6,m = 6) during 
inspiral. However, the higher modes are more sensitive to the order of convergence on the Cauchy 
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side during merger, and hence, the error goes down with a higher convergence order (note that 
we use 8th-order finite difference operators). The phase of (£, m) = (6,6) mode converges at 
sixth-order while the amplitude converges at third order. 

The differences between the waveforms at for resolution h = 0.64M are of order of 0.03% 
in amplitude with a dephasing of 0.002 radians over the course of the evolution. 

It should be noted that the differences in the waves between the two world-tube locations is 
one order of magnitude smaller than the numerical errors in amplitude and phase inherent in the 
Cauchy evolution [|9l . We therefore expect that the dominant error in the waveforms is due to the 
discretisation error during Cauchy evolution. Indeed, using the three resolutions h = 0.96M, h = 
0.80M and h = 0.64M, we find a very similar convergence order for (a) CCE with fixed world- 
tube location, and (b) for a pure Cauchy evolution with the same fixed finite radius extraction 
(compare [9]). 

However, it should also be noted that the (non-convergent) differences between extrapolated 
waveforms and waveforms obtained at J'^ as found in [jH are of the same order as the (convergent) 
discretisation error inherent in the Cauchy evolution for our given resolutions as found in [|60l . 
Further increasing the resolution of the Cauchy evolution will therefore lead to an error that is 
much smaller than the error made by finite radius extrapolation. This highlights the importance 
of the systematic error due to extrapolation from finite radius, which is removed by performing 
computations at J'^ via CCE. 

D. Conservation of energy and angular momentum 

We have used the measured waveforms to compute the initial mass and angular momentum of 
the modeled spacetimes. The initial mass is given by the ADM mass as computed at infinity on 
the initial slice using the methods described in [TTTII . Since the initial slice is conformally flat, the 
initial angular momentum of the spacetime is simply given by (see for example 11761 - 1781 ) 

Jadm = Ci X Pi + Ca X p2 + Si + Sa . (46) 

Here Cj, pi and Sj are the position, the linear momentum and the spin of the i-th black hole. 

After evolution, i.e., after the horizon of the final black hole has equilibriated, we compute the 
mass and spin of the remnant on its apparent horizon using the isolated horizon formalism [r79] - l8TI . 
or alternatively, by using the horizon circumference formula ll82l (the two methods agree within 
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FIG. 7: Convergence of the differences in wave-signal to zero for various test cases. The top left panels 
show the ampUtude and phase differences between world-tube data at rp = lOOM and rp = 200M of 
the (£, m) = (2, 2) mode of the Weyl scalar ^4 for the two resolutions h = 0.96M and h = 0.64M 
where the latter is scaled for second-order convergence. The top right panels show the the same for the 
(£, m) = (4, 4) mode of 1/^4. The lower panels show convergence of the {I, m) = (6, 6) mode. The phase 
error for the high resolution is rescaled for sixth-order convergence while the amplitude error is rescaled for 
third-order convergence. 
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the reported errors as given in Table |l| where we use the values obtained from the circumference 
formula^). In addition, we compute the radiated energy and angular momentum carried by the 
gravitational waves using the results of finite radius computations, from extraction at and 
from an extrapolation of the results at several finite radii. The finite radius computations are 
obtained at an extraction world-tube rr = lOOM , which is typical for current numerical relativity 
simulations (see e.g., [[83l[84ll ). 

The radiated energy and angular momentum flux can be calculated in terms of spin- weighted 
spherical harmonic coefficients of ^4 as derived in [[851 [861 . The radiated energy flux is given by 

2' 



dE , ^ 
— — = lim < > 

dt r^oo IGvr 



dt^^:"" 



^4 



(47) 



and the components of the radiated angular momentum flux read 
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X I £,m+l I X I £,m—l 

Jl>,mVA + Jl-mVi 



dt'ij'r 




f,,m = V£(£+l)-m(m + l). (51) 

Note that we can also compute the expressions for radiated energy and angular momentum in 
terms of the news function by using the fact that in a Bondi-frame, \E'4 is related to the news via 

^4 = 5„iV. (52) 



Hence, we can replace all first time integrals of ^/'4 in ([47[)-(|50|) by A^. 

We find that modes are well resolved only up to (£, m) = (8,8), and hence the sum in 
the expression above is to £ = 8. We have computed the sums for waveforms extracted 
at finite radii and at . The finite radius computations have been performed on extraction 



^ The close agreement between both methods further supports the idea that the final black hole is indeed a Kerr black 
hole. 
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Model 


(ai,a2) = (0,0) 




(01,02) = (0.8,-0.8) 


total ADM mass Madm 


0.99051968 ±2 x 10 


-8 


0.98826246 ± 2 x 10"® 


initial angular momentum Jadm 


0.993300 ± 1.0 X 10- 


16 


0.900312 ± 1.0 X 10-16 


irr. Mass Mf^[„ 


0.884355 ± 2.0 x 10" 


-5 


0.883636 ± 1.0 x 10"^ 


Spin Sf/M^f 


0.686923 ± 1.0 x 10" 


-5 


0.686300 ± 2.7 x 10"^ 


Mass M f 


0.951764 ±2.0 x 10" 


-5 


0.950829 ±3.8 x 10"^^ 


Angular momentum Sf 


0.622252 ± 1.0 x 10" 


-5 


0.620469 ± 1.9 x lO"'' 


Recoil velocity Ukick.ccE 






364.30 ± 0.32 km/s 


Radiated angular momentum Jrad,ccE 


0.370704 ± 6.2 x 10" 


-5 


0.286908 ±2.6 x lO"'' 


Radiated energy E'rad.ccE 


0.038828 ± 2.5 x 10" 


-6 


0.037689 ± 4.3 x 10"^ 


Recoil velocity Ukick,finite,r=ioo 






350.37 ± 1.46 km/s 


Radiated angular momentum Jrad,finite,r=ioo 


0.339725 ±4.0 x 10" 


-4 


0.266654 ±5.7 x lO^^ 


Radiated energy E'rad finite r-ioo 


0.037354 ±3.9 x 10" 


-5 


0.035369 ±3.7 x lO"'' 


Recoil velocity Ukick,extrap 






362.87 ± 1.19 km/s 


Radiated angular momentum Jrad.extrap 


0.370007 ± 6.8 x 10" 


-5 


0.283738 ± 2.1 x lO^^ 


Radiated energy ^^rad.extrap 


0.038715 ±2.2 x 10" 


-6 


0.037189 ± 2.2 x 10~^ 


Residual Vkick,CCE - fkick.extrap 






1.43km/s 


Residual Jrad.CCE — -^rad.cxtrap 


7.0 X lO"'^ 




3.1 X 10-3 


Residual E'rad.CCE — -E'rad.extrap 


1.1 X 10-^ 




5.0 X 10"^ 



TABLE I: Properties of the initial spacetime, of the merger remnant, and of the emitted radiation calculated 
from the appropriate waveform. The reported errors of the extrapolated values refer to the error in the 
extrapolation and the errors of the values at refer to the error from two different extraction world-tubes 
at the same fixed resolution. The error in the initial angular momentum Jadm is at machine precision as it 
is known analytically in terms of the initial parameters of the black holes. All other errors are given by the 
Richardson expanded differences between a medium and a high resolution run as described in [9 |. The last 
three rows display the differences between extrapolated results and results at , which we note is of the 
order of the discretisation error as indicated by the finite radius error-bars. 
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spheres r = {280M, 300M, 400M, 500M, 600M, lOOOM} for the non-spinning binary, and 
r = {lOOM, 250M, 500M} for the spinning binary. Given the radiated energy and angular mo- 
mentum on the three extraction spheres, we do a linear extrapolation by fitting a function of the 
form 

Q(r) =Qo + Qi/r, (53) 

where Q is either i^^ad or Jrad- 

The values for the initial mass Madm» initial angular momentum Jadm^ final (Christodoulou) 
mass Mj, final angular momentum 5/, as well as radiated energy i?rad and radiated angular mo- 
mentum Jrad are reported in Table |I]for both binary systems. The reported errors refer to Richard- 
son expanded differences between a high and a medium resolution run as described in [[91, except 
for the error in the extrapolation and the error in CCE, where for the former we report the error in 
the fit, and for the latter report the differences between two given extraction world-tubes for the 
same fixed resolution. Hence, they do not contain the discretisation error of the evolution, which, 
in all cases, is one order of magnitude larger than the error in the extrapolation and in CCE. 

Table |l] lists the radiated quantities determined from the measured values of ^4. However, 
since we also determine the news function at , we can independently compute i?rad,ccE and 
Jrad.ccE bascd ou in order to cross check the result. We find that the differences in all cases are 
on the order of the reported error bars, hence further validating the results based on CCE. 

Another important point to make is that the differences between extrapolated data and data ob- 
tained at are of the same order of magnitude as the discretisation error in the Cauchy evolution. 
This is indicated by the last three rows of Table |lj The discretisation error of these quantities due 
to the the Cauchy evolution is given by the error-bars of the finite radius quantities. 

Due to conservation of energy and angular momentum, the following relations must hold ex- 
actly 

Sf + Jrad = -^ADM , Mf + Ej-ad = ^ADM • (54) 

The residuals of these relations are reported in Table [II] for the quantities at a finite radius, extrap- 
olated to infinity, and computed at , for both binary systems. 

Note that the residuals in the extrapolated values of the two binaries are of comparable accu- 
racy to the values obtained via CCE. This is surprising, considering the error in the extrapolation 
procedure itself, and demonstrates that extrapolation, even in the most extreme spinning cases, de- 
livers comparable performance. Note, however, that our extrapolation uses large extraction radii 
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Model 



(01,02) = (0,0) 



(01,02) = (0.8,-0.8) 



\Sf + Jrad,finitc - ^ADmI 3.1 X lO'^ 1.3 X lO'^ 

\Sf + Jrad,cxtrap " ^ADm| 1-0 X 10"^ 3 9 ^ 10^3 

\Sf + Jrad,ccE " ^admI 3.4 X 10"^ 7.1 X 10"^ 

\Mf + -Brad.finitc " ^admI 1-4 X 10"^ 2.1 X 10-^ 

\Mf + ^rad.extrap " Madm| 4.1 X 10-^ 2.4 X 10"^ 

\Mf + ^rad.ccE - MadmI 7.2 X 10-5 2.6 x 10"^ 



TABLE II: Residuals in conservation of energy and angular momentum (54i for the non-spinning and 
spinning configuration and for radiated mass E'rad and angular momentum Jrad as calculate at a finite 
radius r = lOOM, extrapolated from finite radii computations to infinity, and by the CCE computation at 
J'^. Both, extrapolated quantities and those directly computed at J'^ satisfy the conservation conditions 
with about the same amount of accuracy. Not using any extrapolation but relying on a pure finite radius 
computation leads to an error that is of about 1-2 orders of magnitude larger than the one found by the other 
two methods. 

compared to common practice in numerical relativity. By using a set of smaller extraction radii 
in the extrapolation, we expect that the errors increase up to one order of magnitude (compare 
ll60ll ). In addition, the systematic error made in the extrapolation will not converge with numerical 
resolution, but only linearly with extraction radius 0{r). This however, puts high computational 
demands on future high accuracy simulations. 

In all cases, we observe that a pure finite radius computation based on an extraction sphere at 
r = lOOM clearly performs worse and has an error that is larger by 1-2 orders of magnitude, 
which, in the case of angular momentum, corresponds to an error of the order of 1.0% — 3.0%. 



E. Recoil velocity 

In addition to energy and angular momentum, gravitational waves also carry linear momentum. 
Depending on the binary configuration, there might be a prefered direction at which radiation, and 
hence linear momentum, is beamed. Due to conservation of linear momentum, this gives rise to a 
gravitational recoil, or "kick" acquired by the merger remnant. An accurate measurement of this 
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is effect is crucial in various fields of astrophysics as this determines whether the final black hole 
is kicked out of its host environment. Clearly, the existence or non-existence of a black hole can 
have a dramatic impact on the host and determines its further evolution. 

Previous studies have already considered this effect in detail, and particularly, the binary con- 
figuration that we consider here belongs to a sequence that has already been examined in ll87l[88l . 
However, these studies estimated the recoil velocity using waveforms extracted at a finite radius. 
In this subsection, we measure this effect where it is defined unambiguously, namely at . 

The radiated linear momentum flux can be calculated in terms of spin-weighted spherical har- 
monic coefficients as derived in [|85l and given by 

dt dt ™\8vrZ.7_ 

f ( 7i,m+l I 7 7^-1, m+1 T 7^+1, m+l\ 1 /rc\ 

X y dt [ai^^i)^ + hi_^ip^ - hi+i^rn+iVi j | > (55) 

dt ™\l67r^y_^ 

X j dt' (^Q,mV^4 + C?Am^4~^'™ + C?£+l,m^4'^^'"') | , (56) 



where 



_ v/(£-m)(£ + m + l) 
<^i,m - , (57) 



1 /(£-2)(£ + 2)(£ + m)(£ + m- 1) 
''"^ ~ (2£-l)(2£ + l) 



= J^y (59) 

_ 1 l ii-2)ii + 2)ii-m)ii + m) 
^^'"^ - £ V (2£-l)(2£+l) • ^^^^ 

Again, we have computed the expression above by including modes up to £ < 8. In addition, 
we calculate this effect using a finite radius computation and a computation at J'^. Given the 
three extraction spheres at r = {lOOM, 250M, 500M} for the spinning configuration, we have 
extrapolated the recoil velocity to infinity by means of Eq. ( [53] ). 

It should be noted that it is crucial to include the appropriate vector integration constant when 
integrating the linear momentum flux in time. As analyzed in [[TOll . there are two important contri- 
butions. The first mainly comes from the initial burst of junk radiation that adds significantly to the 
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measured recoil at the beginning of the simulation. The second is due to the non-zero net-linear 
momentum that the binary had acquired if it had inspiralled from an infinite separation. We have 
measured the correct vector integration constant by applying the method described in [fTOl and cor- 
rected the computation accordingly. We find good agreement between finite radius computations, 
extrapolation to infinity, and computation at . The results are included in Table |lj 

F. Constraint preservation 

We check the preservation of the constraints at the world-tube and at for different radii 
and resolutions. We use the constraint equations as implemented and first tested with linearized 
solutions in [|89l . As a consequence of the Bondi system, there are three constraint equations 

i?oo = 0, i?oi = 0, g^i?oA = 0. (61) 

Analytically, if these constraints are satisfied initially and at the world-tube, they should be satis- 
fied exactly everywhere and at all times. However, in computer simulations, the solution always 
contains numerical error so that the constraints can only be satisfied approximately. Since this er- 
ror is only due to discretisation, it should converge away in the limit of infinite resolution. Indeed, 
this is what we find. 

Fig. [s] shows the convergence of _Roo at and q^Rqa and Rqi at the world-tube F for the non- 
spinning binary. While i?oo appears to be only first-order convergent at , we note that the error 
is already very small and at the order of 10^^° and hence close to machine precision. The same 
applies to the constraint -Rqi which exhibits a convergence exponent between first and second- 
order. The constraint q^Rqa, on the other hand, exhibits an overall second-order convergence 
exponent at the world-tube. The origin of the spurious spikes occuring at the high-resolution run 
between —500 < t < —100 and reducing the convergence order, is currently unkown but may be 
due to the first-order knowledge of U^r at the world-tube. This is subject to future investigation 
and improvements of the algorithm and code. We also note that during merger, the convergence 
of the q^RoA constraint is reduced at the world-tube. This, again, may be due to the first-order 
knowledge of f/^ at the world-tube. However, the convergence of the constraints appeared to be 
problematic already in [|89ll and we emphasize that the differences in the waveforms at J7+ as 
extracted from different world-tube locations converge to zero at approximately second-order and 
independent of the constraints. 
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FIG. 8: The left panel shows the L2-norm of the constraint i?oo at for the two resolutions h = 0.96M 
and h = 0.64M scaled for first-order convergence. The right figure shows the L2-norm of the constraint 
q^RoA at the world-tube T for the two resolutions h = 0.96M and h = 0.64M scaled for second-order 
convergence. The bottom panel shows the same for the constraint i?oi at the world-tube F. All figures refer 
to the non-spinning binary. 

In all cases, the constraints exhibit a large spike at the time when the first non-trivial data, 
particularly the initial burst of junk radiation, reaches the world-tube location. Although the junk 
radiation should satisfy the constraint equations since it arises as part of the solution to the con- 
straint preserving initial data at the initial Cauchy hypersurface, we note that there is principally an 
incompatibility between the conformal flatness assumption of the 3-metric on the initial Cauchy 
hypersurface and the initial data of the 4-metric, J{u = 0) = 0, on the initial null hypersurface. 
This may lead to the observed initial spikes in the constraints. Further investigation is necessary 
to clearly identify the source of error and to improve the construction of characteristic initial data. 
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VI. NOTES ON CHARACTERISTIC CODE PARAMETERS 



There are some subtleties regarding the choice of some code and grid parameters that have been 
found empirically during tests of the extraction method. Here, we collect some of these findings 
for future reference. 

• We use conformally flat initial data on the initial null hypersurface, i.e.J{u = 0) = 0. 
This assumption is in principle incompatible with the proper radiation content of a binary 
black hole spacetime, and it is subject to the same error as on the Cauchy side, where we use 
conformally flat initial data. In practice, the spurious burst of junk radiation at the beginning 
of the simulation quickly leaves the system and does not affect the radiation at the accuracy 
given by the numerics. 

• The radial resolution of the physical coordinate on the characteristic grid should be at least 
as fine at the world-tube as the radial resolution on the Cauchy side. We can estimate the 
compactified radial resolution in terms of the physical resolution via 

which asymptotes to rwt for a; ^ at the world-tube, i.e. 

dr 

a; ^ 0, ^ r„t • (63) 

dx 

It follows that 

Ar|r ~ rwt Ax (64) 

at the world-tube. We can then select an appropriate Ax, so that Ar|r at the world-tube is 
at least as small as the radial resolution Ar on the Cauchy side. 

It turns out that at least ~ 200 radial points are required to resolve the wave properly at 
J^. Since Ai? = 0.64M at the extraction spheres for the highest resolution, the compact- 
ification parameter has to be r^t ~ 256 in order to meet the condition Ar < Ai?. Further 
increasing the characteristic radial resolution while leaving the Cauchy Ai? = const, has no 
noticeable influence on the wave at J^. 

• The linearized approximation for the conformal factor uj is used for the computation of the 
news N. The full non-linear computation results in spurious drifts and does not converge 
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when considering the Schwarzschild test. This is potentially a bug in the existing imple- 
mentation and affects only the news. ^4 is only affected at higher-order which does not 
play a role in practice. However, as the fields at are in the linear regime, the linear 
approximation is a valid simplification. 

• We use the first-order evolution system in the characteristic code, i.e., we use only first eth- 
derivatives. The implementation with second eth-derivatives does not converge and gives 
inconsistent results for the binary waveform when considering different world-tube loca- 
tions. This might also be related to the problem formulated in [[39l l90l . In principle, both 
methods should be equivalent and the non-convergence of the second-order system can be 
attributed to a bug in the code. 

• The news and ^4 are computed using the second-order implementation of the extraction 
module, i.e., using second eth-derivatives. 

• All quantities are decomposed in terms of the complex-valued spin-weighted spherical har- 
monics sY^m- 

VII. CONCLUSIONS 

We have demonstrated results from an implementation of characteristic extraction in numerical 
relativity, applied to the astrophysically relevant test cases of binary black holes. Our implementa- 
tion carries out the characteristic extraction as a post-processing step. As such, the code is general 
purpose and can be applied to boundary data generated by any Cauchy simulation, independent 
of numerical scheme or formulation of the Einstein equations. The code has been tested against 
analytic solutions, and it was found that it exhibits the expected convergence properties, as well as 
gauge-independence, a property lacking in finite radius calculations. Applied to two test cases in- 
volving the inspiral and merger of binary black holes, we have determined the emitted gravitational 
radiation, and found that the CCE method does indeed generate results which are independent of 
the world-tube on which the data has been extracted. We have also determined the residuals in the 
conservation of mass and angular momentum, as well as the recoil velocity (when the black hole 
spins are anti-aligned). We find that the residuals of the CCE method are very similar to those 
obtained by extrapolation to infinity of results obtained using a pure Cauchy evolution, providing 
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an important confirmation of the usefulness of such methods. However, the (non-convergent) dif- 
ferences between both measures is already on the order of the discretisation error of the Cauchy 
evolution, therefore indicating that as the resolution is further increased, the finite radius extrapo- 
lation will be the dominant source of error in the wave-extraction. The error in CCE, as measured 
here for a particular set of resolutions is one order of magnitude smaller, and will converge away 
if the resolution is further increased. 

Characteristic extraction yields results that are free of a number of systematic errors in the mod- 
eUng of the physical problem. Placing the outer boundary at a sufficiently large distance, we are 
able to causally disconnect the artificial boundary from the extraction world-tube. However, there 
is still the possibility of a systematic error in the extent to which there is an unphysical radiation 
content in the initial data. Further, it should be noted that there is, in general, an inconsistency 
between the conformally flat initial data used in the Cauchy computation, and that of J = used 
on the initial characteristic null cone. It is usually assumed that the spurious initial radiation will 
quickly be flushed out of the system in a burst of initial junk radiation. The observed convergence 
of the differences between the radiation signals at different extraction radii, indicates that, after the 
burst of junk radiation, the effect of spurious, unphysical initial data is very small and below the 
discretisation error. 

Future developments will clearly include applying characteristic extraction to a variety of other 
physical systems. It will also be interesting to revisit the problem of Cauchy-characteristic match- 
ing, in which the extraction world-tube and the Cauchy outer boundary coincide and boundary 
information is passed in both directions during the course of the evolution. In this case evolution 
in the Cauchy domain between, say r = lOOM and r = 3000M, is avoided so the method is sig- 
nificantly more efficient; and further one can avoid inconsistencies in the initial data between the 
two domains. Finally, it would be useful to construct initial data that correctly models the physical 
problem; this is, in principle, possible at least for the characteristic domain. 
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